Advances in Applied Acoustics Vol. 1 Iss. 1, November 2012 



Eigenmode Analysis of the Thermoacoustic 
Combustion Instabilities Using a Hybrid 
Technique Based on the Finite Element 
Method and the Transfer Matrix Method 



G. Campa 1 , S.M. Camporeale 2 

Dipartimento di Meccanica, Matematica e Management - 

via Re David 200, 70125 Bari, Italy 

1 g. campa@poliba.it ; 2 camporeale@poliba.it 

Abstract 

Thermoacoustic combustion instabilities affect modern gas 
turbines equipped with lean premixed dry low emission 
combustion systems. In the case of annular combustion 
chambers, experimental test cases carried out on small scale 
test rigs equipped with single burner arrangements fail to 
give adequate indications for the design of the full scale 
combustion chamber, since they are unable to reproduce the 
interaction of the flame fluctuation with the azimuthal 
pressure waves. Therefore there is a large interest in 
developing techniques able to make use of data gathered 
from tests carried out on a single burner for predicting the 
thermoacoustic behavior of the combustion chamber at full 
scale with its actual geometry. 

A hybrid technique based on the use of the finite elements 
method and the transfer matrix method is used to identify 
the frequencies at which thermoacoustic instabilities are 
expected and the growth rate of the pressure oscillations at 
the onset of instability, under the hypothesis of linear 
behavior of the acoustic waves. This approach is able to 
model complex geometries such as annular combustion 
chambers equipped with several burners. Heat release 
fluctuations are modeled through a classical n-x Flame 
Transfer Function (FTF). In order to model the acoustic 
behavior of the burners, the computational domain 
corresponding to each burner is substituted by a 
mathematical function, that is the Burner Transfer Matrix 
(BTM), that relates, one to each other, pressure and velocity 
oscillations at either sides of the burner. Both the FTF and 
the BTM can be obtained from experimental tests or from 
CFD simulations. The use of the transfer matrix permits us 
to take into account parameters, such as the flow velocity 
and the viscous losses, which are not directly included in the 
model. 

This paper describes the introduction of the burner transfer 
matrix in the combustion chamber model. Different 
geometries of combustion chamber and burner are tested. 
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The influence of the parameters characterizing the transfer 
matrix is investigated. Finally the application of the BTM to 
an actual annular combustion chamber is shown. 
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Nomenclature 

Latin: 



A 


cross sectional area 


b 


length 


c 


speed of sound 


F 


Riemann invariant 


G 


Riemann invariant 


i 


imaginary unit 


k 


acoustical wave number 


1 


length 


LHV 


lower heating value 


M 


Mach number 


n 


interaction index 


P 


pressure 


q 


volumetric heat release rate 


Q 


rate of heat release per unit area 


R 


gas constant 


RR 


rate of reaction 


Re 


real part 


s 


thickness 
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t 


time 


T 


temperature 


u 


velocity vector 


X 


position vector 


Greek: 




a 


area ratio 


P 


parameter 


Y 


ratio of specific heats 


5 


Dirac's delta 


C 


pressure loss coefficient 


A 


eigenvalue =-ico 


9 


density 


a 


time lag distribution 


T 


time delay 


CO 


angular frequency 


Subscripts: 




d 


downstream 


eff 


effective 


i 


reference position 


u 


upstream 



Introduction 

The introduction of lean premixed burners in gas 
turbine systems for power generation has determined 
an increase of the risk of thermoacoustic 
instabilities [1,2,3]. This phenomenon is caused by self- 
sustained pressure oscillations, which arise from 
coupling of heat release fluctuations and acoustic 
waves propagating within the combustion system. 
These instabilities are dangerous because can lead to 
strong vibrations of the system with consequent 
structural damages. The phenomenon is known since 
a long time and several models and techniques have 
been proposed over the years to have a better 
comprehension of its causes and to predict which 
operating conditions are more prone to cause it. 

Quite often low-order models [1,4,5,6,7] are proposed 
to represent the combustion instability in the 
combustion system modeled as a network of acoustic 
elements (and so the as much used name of Acoustic 
Networks), where each element corresponds to a 
component, such as a duct, a nozzle, a burner, etc... 



Generally, each component represents a one- 
dimensional propagation of the acoustic waves, 
mainly within the limits of the linear acoustic and 
harmonic time dependence (icot). In case of the 
application to annular chambers, two-dimensional 
components are also proposed to model both axial and 
azimuthal waves [8,9]. These tools may be used to 
study the pressure oscillations by dividing the 
combustor into several zones and representing the 
acoustic field in each zone. Hence, a combustor system 
is defined as a series of subsystems, using 
mathematical transfer function matrices to connect 
these lumped acoustic elements one to each other, so 
providing the continuity of acoustic velocity and 
pressure across each zone. The unknowns are acoustic 
pressure and velocity fluctuations (p' and u 
respectively) at the ports of each network element. In 
this scheme flame is concentrated in one of the lumped 
elements and located at the beginning of the 
combustion chamber. Flame transfer function (FTF) is 
written as the ratio of heat release rate oscillation to 
inlet velocity fluctuation, where both of them are 
normalized by their corresponding time-averaged 
values. These models are very helpful to understand 
the instability mechanisms, but when complex three- 
dimensional geometries are involved, these models 
become unable to model the propagation of the 
acoustic waves adequately. 

LES (Large Eddy Simulation) codes are proposed by 
several authors [10,11,12,13,14] to investigate the 
phenomenon of combustion instability and matching 
pressure oscillations with turbulent combustion 
phenomena, considering the dynamically relevant 
scales during the resolution, related to fluid dynamic 
phenomena and the kinetic chemistry. The 
disadvantages are concerned with the enormous 
computational efforts required for the resolution. 

A finite element method (FEM) approach may also be 
used to solve the acoustic problem in three- 
dimensional geometries [15,16,17]. Although some 
authors proposed the acoustic analysis in the time 
domain in order to model the onset of the 
perturbations and how they increase reaching the limit 
cycle [15], in [16,17] the complex eigenfrequencies of 
the system are detected to ascertain if the 
corresponding modes are stable or unstable [17]. This 
approach numerically solves the differential equation 
problem converted in a complex eigenvalue problem 
in the frequency domain. The eigenvalue problem is 
nonlinear and it is solved by means of a linearization 
under the hypothesis of small oscillations. 
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In this paper a hybrid technique is described, since the 
combination of lumped elements and three- 
dimensional elements is achieved. The numerical 
procedure, able to detect the complex eigenvalues of 
the analyzed system in presence of a flame transfer 
function and a transfer matrix function, is 
implemented within commercial software, COMSOL 
Multiphysics, based on the finite element method. The 
main idea is to introduce in the FEM framework the 
transfer matrices of the low-order models to model the 
burners. Each burner is described not as a 
computational domain, but as a mathematical function, 
which can be obtained analytically or experimentally. 
Additionally, this idea permits us to take into account 
the influence of parameters, such as the mean flow 
and the viscous losses, which are not considered in the 
FEM framework. A strong reduction in the 
computational time is also obtained. 

The goal of this study is to provide adequate 
indications in the design stage of a full scale 
combustion chamber. In fact, there is a great interest 
on developing methods able to make use of the data 
gathered from experimental or numerical tests carried 
out on a single burner for predicting the thermo 
acoustic behavior of a practical combustion chamber. 

In the first section the mathematical model is 
introduced. In the second section the application of the 
transfer matrices is described and the applications to 
different geometries of the combustion chamber are 
shown: cylindrical and annular configurations. The 
influence of parameters characterizing the flame 
model and the transfer matrix is also investigated. In 
the third section the transfer matrix method is applied 
to a practical annular combustion chamber. 

Mathematical Model 

As regards gas turbine combustion systems, the flow 
velocity is generally far below the sound velocity, both 
in the plenum placed at the exit of the compressor and 
in the combustion chamber, while the flow velocity is 
not negligible in some limited areas, such as conduits 
of single burner or multiple burners that connect the 
plenum to the chamber. These areas, in terms of 
propagation of pressure waves, can be treated as 
"compact" elements that can be modeled by means of 
specific transfer function matrices, obtained 
experimentally or numerically through CFD or 
aeroacoustics codes [18,19,20,21]. The hybrid method 
proposed here makes a combined use of 



■ FEM for the analysis of the propagation of the 
acoustic wave in three dimensional domains 
with negligible flow velocity, 

■ Transfer matrices to model the acoustic waves 
in the ducts of the burners, where there is a 
one-dimensional propagation of the waves, in 
presence of not negligible levels of flow 
velocity. 

In this section the mathematical model adopted for the 
FEM analysis of the volumes with negligible flow 
velocity, even in presence of heat release fluctuations 
and temperature gradients, is described. 

In the acoustic analysis each variable is assumed to be 
composed of a steady mean value (identified by an 
overbar) and a small perturbation (identified by a 
prime): p(x,t) = p(x) + p'(x,t) and u (x,f) = u(x) + u'(x,f) . 
The flow velocity is considered negligible in 
comparison with the sound velocity, within the 
computational domain, so u(x) =0 in the whole 
domain. The effects of viscous losses and thermal 
diffusivity can also be neglected, while the 
temperature can vary within the domain under the 
simplifying hypothesis that the fluid can be 
considered as an ideal gas, which means that specific 
heats are supposed constant and the variation of the 
composition caused by combustion does not produce a 
variation of the number of moles. From these 
hypotheses, in presence of heat fluctuations, the 
linearized equations for the perturbations can be 
obtained [8,22], giving the inhomogeneous wave 
equation 



1 d 2 p' 
c 2 dt 2 



-pv 



4vy 



y-ldq' 
c 2 dt ' 



(1) 



where c is the speed of sound and q' is the fluctuation 
of the heat input per unit volume. The term at the RHS 
of Eq.(l) shows that the rate of non-stationary heat 
release creates a monopole source of acoustic pressure 
disturbance. Considering that mean flow velocity is 
neglected, no entropy wave is generated and the 
pressure fluctuations are related to the velocity 
fluctuations by 



dt p 



(2) 



This work is conducted in the frequency domain, so 
that each variable fluctuation is expressed as a 
complex function of time: 
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p'= Re( p(x) exp(icot)) 



(3) 



where the hat refers to a complex variable and the 
same is for the other variables. Introducing this 
condition into Eq.(l), the Helmholtz equation can be 
obtained: 



X 2 , 
c 



pV- 



P 



c 



(4) 



where X — —ico is the eigenvalue, being co a complex 
variable, comprising a real part that gives the 
frequency of oscillations and an imaginary part that 
gives the growth rate at which the amplitude of 
oscillations increases per cycle. Eq.(4) shows a 
quadratic eigenvalue problem which can be solved by 
means of an iterative linearization procedure. 
COMSOL Multiphysics uses the Arpack Fortran as 
numerical routines for large-scale eigenvalue 
problems. It is based on a variant of the Arnoldi 
algorithm, called the implicit restarted Arnoldi 
method [23]. An iterative procedure based on a 
quadratic approximation around an eigenvalue 
linearization point X 0 is adopted. Such a procedure is 
speeded up by using, as approximate starting 
eigenvalue, the value obtained through the analysis of 
the system without heat release fluctuations. The 
solver reformulates the quadratic eigenvalue problem 
as a linear eigenvalue problem of the conventional 
form Ax = XBx , and, iteratively, updates the 
linearization point until the convergence is reached. 
The software approximates the nonlinear model with a 
linear model, and the approximation is valid only 
when the solution is close to the linearization point. 
This means that not always the obtained solution is the 
exact solution, but it could be a spurious solution, so 
that a manual re-initialization is needed. 

Transfer Matrix Method 

The transfer matrix is a 2x2 matrix whose coefficients 
can relate the fluctuations of acoustic pressure p' and 
velocity u' at one junction, or port of the element, to 
the fluctuations of acoustic pressure and velocity on 
the other junction of the element. Instead of the 
physical acoustic variables, the Riemann invariants F 
and G are usually adopted. Riemann invariants are 
related to acoustic pressure and velocity by 



G = 



K pc 



(5) 



and can be thought as waves propagating, respectively, 
in the downstream and upstream direction. For a 
simple duct of uniform cross section, under the 



hypotheses of plane waves without mean flow and 
dissipative effects, the transfer matrix is given by [24]: 



f Z7 A 



\ G dJ 



r exp(ikl) 0 
0 exp(-ikl) 



v 



f Z7 \ 



\ G uJ 



(6) 



where k = col c is the wave number and I is the length 
of the duct. 

Test Cases 

At the beginning, some preliminary application tests 
are carried out on a duct with uniform cross-section 
area, uniform temperature and no mean flow. This 
scheme is characterized by an open-end inlet (p'=0) 
and a closed-end outlet (u'=0), see Fig. 1. A planar 
flame sheet is supposed to be located at the abscissa 
x=b. Such a scheme is the same proposed by Dowling 
and Stow [8]. In Fig. 2, the scheme A represents the 
computational domain of conventional "one-piece" 
duct, while the computational domain of the scheme B 
is divided in two pieces connected by means of the 
transfer matrix, given in Eq.(6), where / represents the 
length of the piece of the duct substituted by the 
transfer matrix. In this one-dimensional case, the 
planar wave hypothesis can be adopted and only the 
abscissa x along the duct is used to describe the 
variation of acoustic pressure and velocity in the duct. 



plenum 
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FIG. 1 SCHEME OF A STRAIGHT DUCT WITH UNIFORM CROSS 
SECTION 




FIG. 1 COMPUTATIONAL GRIDS OF THE STRAIGHT DUCTS 
WITH UNIFORM CROSS SECTION. DUCT A IS WITHOUT BTM, 
DUCT B IS WITH BTM 
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Test Case 1 

In the first test case the transfer matrix method is 
applied to the described scheme without flame 
fluctuations. The results obtained by means of the 
application of the proposed method are compared 
with those obtained from the same scheme without the 
interruption under the same operating conditions. 



Q\t) - 



Ml 
r -i 



(8) 




FIG. 2 MODESHAPE CORRESPONDING TO THE FIRST 
EIGENMODE OF THE STRAIGHT DUCT WITH UNIFORM CROSS 
SECTION (FIG. 2) WITH UNIFORM TEMPERATURE. RED LINE 
FOR DUCT B, BLUE LINE FOR DUCT A 

For comparison, in Fig. 3, the mode shapes 
corresponding to the first eigenfrequency are shown 
for both the schemes, A and B. It appears that there is 
a very good agreement between the results obtained 
with transfer matrix and those obtained from the 
original "one-piece" duct. The correspondence of the 
data is highlighted in the zoom of acoustic pressure 
path. The first eigenfrequency of this system results 
equal to 21.4 Hz in both cases. Pressure and velocity 
fluctuations on one side are well linked to the 
corresponding on the other side of the interruption. 
The results confirm that, in such a simple geometry, 
without changes of cross section geometry and area, 
where only plane waves are present and no tri- 
dimensional effects, the transfer matrix does not 
determine any alteration of the modeshapes. 

Test Case 2 

In the second test case the transfer matrix method is 
applied to the described scheme including flame 
fluctuations. Following the approach used in [8], the 
rate of heat release fluctuations per unit volume q'(x,t) 
is assumed to be related to the local velocity u' 
upstream the flame zone with a time delay t, by: 



where Q'(t) is the rate of heat input per unit area of the 
cross section of the duct and subscript 1 denotes 
conditions just upstream of this region of heat input: 
«]'(/) = u'(x[,t), where x x -s-b/2 (see Fig. 1). Eq.(7) 
relates the fluctuations of heat input rate per unit 
volume, Cj'(x,t), to the fluctuations of heat input rate 
per unit area of the cross section, Q'(t), through the 
Dirac's delta d(x-b). The non dimensional parameter f> 
gives a measure of the coupling between heat release 
and velocity fluctuations. In the FEM eigenvalue 
analysis, heat release fluctuations are assumed to 
occur in a thin volume with thickness s and the Dirac's 
delta 5(x-b) can be roughly approximated by 



S(x-b)-. 



0 x<b-s/2 
l/s b-sl2<x<b + sl2 
0 x>b + sl2 



(9) 



Considering Eq.(3), after some algebraic steps, Eq.(4) 
becomes 

-—rk 1 p - — V ■ (Vp) = [- pS{x - b)(-A)u(b') exp(Ar)] (10) 
pc p 

that is the governing equation to be solved in the FEM 
code. 

Fig. 4 shows the modeshapes for the first axial mode, 
considering that blue line is related to the whole "one 
piece" system, whereas red line is related to the 
system with transfer matrix. 




q\x,t) = Q\t)S{x-b) 



(7) 



FIG. 3 MODESHAPES CORRESPONDING TO THE FIRST 
EIGENMODE OF THE DUCT IN FIG. 2 WITH Ti/Ti=2.14 AND 
B/L=0.3. RED LINE FOR DUCT B, BLUE LINE FOR DUCT A. 

In Fig. 4 the position of the flame can be recognized by 
the presence of a knee in the modeshape. Modeshapes 
with and without transfer matrix are nearly coincident 
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and the differences can be considered negligible. Table 
1 highlights the good agreement of the results with 
and without the transfer matrix. 

TABLE 1 VALUES OF THE FIRST EIGENFREQUENCIES 
CORRESPONDING TO THE EIGENMODES IN FIG. 4. 





b/L = 0.3 


b/L = 0.5 


with transfer matrix 


(34.3 + 1.3i) Hz 


(31.6 + l.li) Hz 


without transfer 
matrix 


(34.8 + 1.3i)Hz 


(32.0 + l.li) Hz 



Test Case 3 

A little more realistic quasi one-dimensional 
combustor is now presented. The analyzed scheme is 
characterized by three cylindrical ducts, representing 
plenum, burner and combustion chamber, respectively. 
In this case, the burner has a diameter lower than that 
of plenum and combustion chamber, which, for 
simplicity, has the same diameter (Fig. 5), with a ratio 
between the length of the burner, y, and the overall 
length, L, is equal to 0.1. The boundary conditions at 
the inlet and at the outlet are assumed as acoustically 
closed ends, u '(0) =0 and u '(L) =0. 

Fig. 6 shows the criteria adopted to substitute the 
burner with a transfer matrix: scheme A represents the 
original geometry, while schemes B and C represent 
two alternative geometries where a transfer matrix is 
used to model the removed volume of the burner. The 
difference between the schemes B and C concerns the 
surfaces adopted as interfaces between the transfer 
matrix (ID model) and the domain discretized by FEM 
(3D model). Such surfaces are colored in red in Fig. 6. 
In the scheme C both the interface surfaces are circles 
with the same diameter of the burner, whereas in the 
scheme B the interfaces have the same diameter of 
plenum and combustion chamber. 

Heat release fluctuations are modeled by Eq.(10). 



Heat release zone 



premixer 



plenum 



|12 

T 2 Combustion chamber 
c 3 



X — h 

L 



FIG. 4 SCHEME OF A STRAIGHT DUCT WITH VARIATION OF 
SECTION. 

Fig. 7 and Table 2 show that the application of 



analytical transfer matrix to the system with variation 
of section yields very good results both in modeshapes 
and in identifying the complex eigenfrequencies. 



scheme f ! 




scheme A 



FIG. 5 DIFFERENT COMPUTATIONAL DOMAINS AND 
INTERFACES TO THE BTM (PARTICULAR). SCHEME A 
REPRESENTS THE ORIGINAL GEOMETRY, SCHEMES B AND C 
WITH BTM. 

TABLE 2 VALUES OF THE EIGENFREQUENCIES. 



Scheme 


1 st mode 


2 nd mode 


A 


(77.4 + 4.2i) Hz 


(212.8 - 0.3i) Hz 


B 


(79.7 + 3.7i) Hz 


(213.4 - 0.4i) Hz 


C 


(77.2 + 4.2i) Hz 


(212.7 - 0.3i) Hz 



Scheme C gives results better than B, because the 
effects of the radial propagation of the waves 
occurring in the original scheme A are well 
reproduced, as it is shown in Fig. 8, which clearly 
shows that scheme C is able to capture the radial 
component of the acoustic waves in proximity of the 
burner, whereas scheme B does not recognize such 
effects. 
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From this test case, it appears that if interfaces are 
correctly identified, the model with transfer matrix is 
able to reproduce very accurately the original system. 
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FIG. 6 MODESHAPES CORRESPONDING TO THE FIRST (A) AND 
SECOND (B) EIGENMODE OF THE DUCT WITH VARIATION OF 
SECTION (FIG. 6) AND T2/Tj=2. RED LINE FOR SCHEME C, 
GREEN LINE FOR SCHEME B, BLUE LINE FOR SCHEME A 
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FIG. 7 ACOUSTIC PRESSURE FIELD FOR SCHEMES A, B AND C 

Test Case 4 

In the following test the burner is modeled by means 
of the transfer matrix proposed both by Fanaca [20] 
and Alemela [21] in order to model the burner used 
during their experiments at Technische Universitet of 
Munich. The transfer matrix is obtained assuming one- 
dimensional flow with low Mach number within a 
"compact element" represented by a duct having short 
length (compared to the wavelength), variable cross 
section and pressure losses. Linearizing the mass and 
momentum conservation equations, the following 
equations can be obtained 



-M +u 



ICO 



"Jeff 



pc 



Mu 



P 
pc 



(11) 



-£M d u d =0 (12) 



where M is the Mach number. In Eq.(12) the effective 
length l-ff given by 



/ = f Au dx 
eS U(x) ' 



(13) 



takes into account the inertia of the air mass in the 
duct, assuming that the Mach number is sufficiently 
low and that the effect of air compressibility can be 
neglected. The coefficient C gives the acoustical 
pressure losses and is generally close to the mean flow 
pressure loss coefficient, due to friction and flow 
separation g = 2Ap/(pu 2 ) . Using effective length L/f 
and pressure loss coefficient Q, neglecting higher order 
Mach number terms, the transfer matrix of a compact 
element is obtained from Eq.(ll) and Eq.(12) 



P_ 
pc 



1 M u -aM d (l + £)-ikl^ 
aM u -M d a + M d ikl eff 



JL 

pc 



(14) 



where a = A u l A d is the area ratio and k = col c is the 
wave number. 

The transfer matrix of a burner modeled as a compact 
element by Eq.(14), is applied to the same one- 
dimensional combustion system examined in Test 
Case 3 (Fig. 5) with addition of pressure losses and not 
negligible Mach number. Fig. 9 shows the influence of 
the Mach number on the first eigenmode. In this way 
it is possible to carry on this procedure both in design 
stage and in check stage. In the design stage, one can 
simulate the process by varying the operating 
parameters in order to obtain the desired specimens. 
In the check stage, one can verify that the solutions 
adopted in order to avoid thermoacoustic combustion 
instabilities succeed in their aim, introducing the 
obtained experimental results or modeling these 
solutions inside the finite element code. 

Additionally, the influence of some parameters (Uff and 
C), characterizing the transfer matrix Eq.(14), on 
frequency and growth rate of the modes is conducted. 




FIG. 8 MODESHAPES CORRESPONDING TO THE FIRST 
EIGENMODE OF THE STRAIGHT DUCT WITH VARIATION OF 
SECTION SHOWN IN FIG. 5, AND T2/Ti=2, Z=1.2. 
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FIG. 9 INFLUENCE OF LOSS COEFFICIENT AND EFFECTIVE 
LENGTH ON NORMALIZED FREQUENCY (A) AND GROWTH 
RATE (B) OF THE FIRST EIGENMODE 

Fig. 10a shows that frequency is influenced only by the 
effective length of the burner. In fact when the length 
of the burner increases, the frequency of the system 
decreases, but it is not influenced by the variation of 
the loss coefficient C,. On the other hand, Fig. 10b 
shows that growth rate is influenced both by the loss 
coefficient C and by the effective length. The mode is 
stable for this configuration, but it becomes more 
stable when C increases, determining a larger damping. 
These results show that it is possible to avoid 
combustion instabilities correctly designing the shape 
of the burners taking into account the geometrical 
constraints of the combustion chamber. In this way a 
correct setting of the fluid dynamic conditions inside 
the burners have to be well considered, in order to 
avoid important decrement of static pressure. 

Test Case 5 

In this test case the configuration is the same of the 



previous test, except the shape of the burner, which is 
conical (Fig. 11), with the ratio of cross-section area of 
the entrance to the burner and of the plenum equal to 
0.562, the ratio of cross-section area of the exit of the 
burner and of the combustion chamber equal to 0.125. 
The burner is substituted by the transfer matrix of 
Eq.(14) and heat release fluctuations are considered. 
Two different conditions are analyzed to observe the 
influence of the Mach number and the influence of the 
pressure loss coefficient, C 
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FIG. 10 SCHEME OF A STRAIGHT DUCT WITH VARIATION OF 
SECTION AND CONICAL BURNER. 
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FIG. 11 VARIATION WITH MACH NUMBER OF FREQUENCY 
AND GROWTH RATE OF THE FIRST EIGENMODE FOR THE 
STRAIGHT DUCT WITH CONICAL BURNER (FIG. 11). 

Fig. 12 shows the influence of the Mach number on 
frequency and growth rate with a fixed value of C As 
in the previous test, growth rate decreases when Mach 
number increases. On the other hand, frequency 
increases with Mach number. Fig. 13 shows the 
influence of the pressure loss coefficient C with a fixed 
value of the Mach number. The trend is similar to the 
one obtained varying the Mach number. 
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X 5 = 0 




0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 

FIG. 12 VARIATION OF FREQUENCY AND GROWTH RATE OF 
THE FIRST EIGENMODE FOR Z =0 TO 3 WITH AN INCREMENT 
OF 0.5 FOR THE STRAIGHT DUCT WITH CONICAL BURNER 
(FIG. 11). 

Test Case 6 

This test case concerns with a more complex geometry, 
that is the annular configuration examined by 
Pankiewitz and Sattelmayer [15]. The combustor is 
characterized by a diffusion chamber ring (plenum) 
and an annular combustion chamber connected by 12 
cylindrical burners, Fig. 14. The mean diameter is 
0.437 m, the external diameter of the plenum is 0.540 m 
and of the combustion chamber is 0.480 m. The length 
of the plenum is 0.200 m and of the combustion 
chamber is 0.300 m. Each burner has a diameter of 
0.026 m and a length of 0.030 m. Temperature in the 
combustion chamber is 2.89 times the temperature on 
the plenum. Flame is assumed to be concentrated in a 
narrow zone at the entrance of the combustion 
chamber, as shown in Fig. 14. Closed-end inlet and 
outlet are assumed as boundary conditions (u'=0). 



plenum premixer 



Ci 


L- 




< 



Combustion 
chamber 




FIG. 13 GEOMETRY OF THE ANNULAR COMBUSTION 
CHAMBER. 

In this case, following the method proposed by Waltz 
et al. [25], and by Forte et al. [26], only a quarter of the 
whole annular system, divided by two meridian 
planes as shown in Fig. 15, is considered as 
computational domain and discretized in the FEM 



analysis in order to reduce the computational effort. 
The analysis of one quarter of the whole geometry, 
requires two steps: in the first one a closed end (u'=0) 
is assumed as boundary condition on both the 
meridian planes (symmetrical boundary conditions); in 
the second step a closed end is assumed on one 
meridian plane and an open end (p'=0) at the other 
plane (asymmetrical boundary conditions). The 
eigenmodes of the original system are obtained as the 
ensemble of the eigenmodes obtained in the first and 
in the second step. Through this approach a smaller 
domain can be simulated, able to provide the same 
eigenvalues and eigenmodes of the original geometry. 
In this test the burners are modeled by means of 
analytical transfer matrices. 

The obtained results for the first four modes are 
provided in Table 3. Eigenmodes are denoted with the 
nomenclature (l,m,n), where /, m and n are, 
respectively, the orders of the pure axial, 
circumferential and radial modes. The corresponding 
eigenmodes are shown in Fig. 15. 

TABLE 3 VALUES OF THE NORMALIZED FREQUENCIES. 



Mode Shape 


(1,0,0) 


(0,1,0) 


(1,1,0) 


(0,2,0) 


Numerical 
Results 


0.242 


0.350 


0.576 


0.658 




(1,1,0) 



FIG. 14 FIRST FOUR ACOUSTIC EIGENMODES OF THE 
COMBUSTOR. 

Fig. 16 clearly shows that the use of the transfer matrix 
does not alter the acoustic pressure field around the 
burners. 




FIG. 15 ACOUSTIC PRESSURE FIELD WITH (LEFT) AND 
WITHOUT (RIGHT) TRANSFER MATRIX 
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The influence of the pressure loss coefficient is 
examined by varying the time delay % for constant 
Mach number and f>. Two different values of C are 
considered: C = 0 and C, = 2 and the paths of normalized 
frequency and growth rate are detected for the first 
circumferential mode in the plenum (Fig. 17) and the 
first circumferential mode in the combustion chamber 
(Fig. 18). 

Fig. 17 and Fig. 18 show the great importance of the 
pressure loss coefficient and, generally, of the 
damping effect which can be applied. An increase of £ 
determines a reduction in the growth rate and so an 
increase in the stability condition. The effect on the 
frequency has not a great importance, because its 
variation is negligible. 
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(b) 



FIG. 16 VARIATION WITH THE TIME DELAY OF NORMALIZED 
FREQUENCY (A) AND GROWTH RATE (B) OF THE FIRST 
CIRCUMFERENTIAL MODE IN THE PLENUM (FIG. 15) 

In the light of these results, it is worth noting the 
combined effect of time delay and pressure loss 
coefficient. A combination of these two elements, 
together with the effect of the Mach number, can be 
helpful to design the burner with parameters able to 
lead to stable conditions inside the combustion 
chamber. 



Application to a 
Combustion Chamber 



Practical Annular 



The idea to model the burners with the transfer matrix 
method is extended to a practical application. The 
annular combustion chamber by Ansaldo Energia is 
examined. An annular plenum is placed at the exit of 
the axial compressor and is connected to the annular 
combustion chamber through 24 burners, whose 
details can be found in [27]. Fig. 19 shows the 
geometrical configuration and the boundary 
conditions. 
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FIG. 17: VARIATION WITH THE TIME DELAY OF 
NORMALIZED FREQUENCY (A) AND GROWTH RATE (B) OF 
THE FIRST CIRCUMFERENTIAL MODE IN THE COMBUSTION 
CHAMBER (FIG. 15). 
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FIG. 18 COMPUTATIONAL DOMAIN AND BOUNDARY 
CONDITIONS OF THE ANNULAR COMBUSTION CHAMBER BY 
ANSALDO ENERGIA 
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In the acoustic analysis, the volume corresponding to 
the duct of each burner is removed and substituted by 
a transfer matrix, so that the upstream port of each 
matrix is the exit from the plenum and the 
downstream is the inlet of the combustion chamber. 
Operating conditions are taken from experimental 
data and from steady flow Reynolds Averaged Navier 
Stokes (RANS) simulations. A spatial distribution of 
the temperature is defined using the results of the 
RANS simulations, which have been imported into the 
3D FEM solver. The heat release fluctuations are 
modeled using the flame model 

gx) = _ wM /(;-r(x)) (15) 
q(x) u t 

In this configuration the reference position i for the 
acoustic velocity corresponds to the combustion 
chamber inlet, which is also the burner transfer matrix 
interface. The mean value of the heat release q is 
expressed as 

q = RR(x) ■ LHV (16) 

where RR represents the spatial distribution of Rate of 
Reaction [kmol/m 3 s], LHV is the Lower Heating Value 
[J/kmol] of the fuel. The volumetric heat release model 
defined in Eq.(16) is transferred into the governing 
equation. In doing so heat release law is applied inside 
the whole combustion chamber domain, and not in a 
specific thin domain at the beginning of the 
combustion chamber as done in the previous models. 
All the information concerning this modeling in the 
FEM code are shown in a previous work of ours [28]. 

An investigation concerning the influence of the 
pressure gradient Ap inside the burners is carried out. 
The influence of Ap is considered inside the burner 
transfer matrix, Eq.(14), and particularly through the 
pressure loss coefficient C,. The linear flame response, 
Eq.(15), with the spatial distribution of the Rate of 
Reaction is used and the complex eigenfrequencies are 
detected and shown in Fig. 20. The time delay is 
assumed to be constant, t = 7 ms. Three different 
values of £ are considered: the pressure loss coefficient 
corresponding to the design burner (C = 0.98), the 
pressure loss coefficient corresponding to a larger 
pressure gradient (C = 2.27), and a zero pressure loss 
coefficient. 




. 20 t i i i i i i L 



0.5 1 1.5 2 2.5 3 3.5 4 4.5 

Normalized Frequency 

FIG. 19 COMBUSTION CHAMBER MODES FOR DIFFERENT 
VALUES OF Z WITH THE SPATIALLY DISTRIBUTED FLAME, 
EQ.(15). T=7MS. 

The results in Fig. 20 show damping effects produced 
by the pressure losses: the increase of C determines a 
general decrease in the growth rate, with some modes 
moving from the unstable area to the stable one. Not 
all the modes are subjected to a change in the 
frequency and in the growth rate, some of them are 
not influenced. Modes 3, 6 and 7 have the largest 
variation in the growth rate, moving from an unstable 
condition without damping (C = 0) to a stable one 
(with C = 0.98 and £ = 2.27). The causes of these results 
can be found in the fact that Modes 3 and 7 are axial 
modes that involve propagation of the pressure waves 
through the burners, whereas Mode 6 is an azimuthal 
one involving propagation of the pressure waves over 
the entire combustor system. In general, the modes for 
which there is an interaction between the waves in the 
plenum and the waves in the combustion chamber are 
characterized by a stronger influence from the 
pressure loss coefficient C- On the other hand, C has no 
influence on the modes for which the wave interaction 
is negligible, such as for Mode 4, and on the modes for 
which the fluctuations are negligible inside the 
combustion chamber, such as for Mode 8. 

Conclusions 

A hybrid technique, based on the combined use of a 
finite element method and acoustic transfer matrices, 
is proposed for the simulation of the combustion 
instability. In order to establish proof of concept, the 
method is applied to several test cases, from simple 
cylindrical ones to an annular one and finally to an 
industrial annular configuration. The substitution of a 
piece of geometry can be useful when this element has 
a complex shape, so that it is possible to reduce the 
computational efforts. The search of the instability 
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conditions by means of this method yields good 
results. The modes are caught very well, and the 
possible differences are due to tridimensional effects 
which at high frequencies become more and more 
important despite to the plane wave theory, that can 
be correctly applied to lower frequencies. 

This technique is appropriate to treat complex 
geometries based on real dimensions and shapes, 
which are really difficult to be treated by means of 
analytical methods or acoustic network methods. In 
this analysis experimental transfer matrix functions 
and flame transfer functions are successfully applied. 

In conclusion, this hybrid technique appears to be a 
simple tool for the analysis of thermoacoustic 
combustion instability both in design stage and in 
check stage, taking the hypothesis of linear acoustic 
waves. The application of this method in the design 
stage permits to keep an useful instruments to detect 
the best operating parameters and the best geometrical 
shapes and dimensions of the examined system, 
without the necessity of facing long and expensive 
experimental tests. In the check stage this method 
permits to verify the correctness of the assumptions 
made in the design stage detecting the stable and 
unstable conditions comparing them to the effective 
operating conditions. 
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